
/* fix_fft.c - Fixed-point in-place Fast Fourier Transform */

/*

All data are fixed-point short integers, in which -32768

to +32768 represent -1.0 to +1.0 respectively. Integer

arithmetic is used for speed, instead of the more natural

floating-point.

For the forward FFT (time -> freq), fixed scaling is

performed to prevent arithmetic overflow, and to map a 0dB

sine/cosine wave (i.e. amplitude = 32767) to two -6dB freq

coefficients. The return value is always 0.

For the inverse FFT (freq -> time), fixed scaling cannot be

done, as two 0dB coefficients would sum to a peak amplitude

of 64K, overflowing the 32k range of the fixed-point integers.

Thus, the fix_fft() routine performs variable scaling, and

returns a value which is the number of bits LEFT by which

the output must be shifted to get the actual amplitude

(i.e. if fix_fft() returns 3, each value of fr[] and fi[]

must be multiplied by 8 (2**3) for proper scaling.

Clearly, this cannot be done within fixed-point short

integers. In practice, if the result is to be used as a

filter, the scale_shift can usually be ignored, as the

result will be approximately correctly normalized as is.

Written by: Tom Roberts 11/8/89

Made portable: Malcolm Slaney 12/15/94 malcolm@interval.com

Enhanced: Dimitrios P. Bouras 14 Jun 2006 dbouras@ieee.org

*/


/*

Henceforth "short" implies 16-bit word. If this is not

the case in your architecture, please replace "short"

with a type definition which *is* a 16-bit word.

*/

/*

Since we only use 3/4 of N_WAVE, we define only

this many samples, in order to conserve data space.

*/

/*
const myint2  Sinewave[N_WAVE-N_WAVE/4] =
{

    0, 201, 402, 603, 804, 1005, 1206, 1406,

    1607, 1808, 2009, 2209, 2410, 2610, 2811, 3011,

    3211, 3411, 3611, 3811, 4011, 4210, 4409, 4608,

    4807, 5006, 5205, 5403, 5601, 5799, 5997, 6195,

    6392, 6589, 6786, 6982, 7179, 7375, 7571, 7766,

    7961, 8156, 8351, 8545, 8739, 8932, 9126, 9319,

    9511, 9703, 9895, 10087, 10278, 10469, 10659, 10849,

    11038, 11227, 11416, 11604, 11792, 11980, 12166, 12353,

    12539, 12724, 12909, 13094, 13278, 13462, 13645, 13827,

    14009, 14191, 14372, 14552, 14732, 14911, 15090, 15268,

    15446, 15623, 15799, 15975, 16150, 16325, 16499, 16672,

    16845, 17017, 17189, 17360, 17530, 17699, 17868, 18036,

    18204, 18371, 18537, 18702, 18867, 19031, 19194, 19357,

    19519, 19680, 19840, 20000, 20159, 20317, 20474, 20631,

    20787, 20942, 21096, 21249, 21402, 21554, 21705, 21855,

    22004, 22153, 22301, 22448, 22594, 22739, 22883, 23027,

    23169, 23311, 23452, 23592, 23731, 23869, 24006, 24143,

    24278, 24413, 24546, 24679, 24811, 24942, 25072, 25201,

    25329, 25456, 25582, 25707, 25831, 25954, 26077, 26198,

    26318, 26437, 26556, 26673, 26789, 26905, 27019, 27132,

    27244, 27355, 27466, 27575, 27683, 27790, 27896, 28001,

    28105, 28208, 28309, 28410, 28510, 28608, 28706, 28802,

    28897, 28992, 29085, 29177, 29268, 29358, 29446, 29534,

    29621, 29706, 29790, 29873, 29955, 30036, 30116, 30195,

    30272, 30349, 30424, 30498, 30571, 30643, 30713, 30783,

    30851, 30918, 30984, 31049, 31113, 31175, 31236, 31297,

    31356, 31413, 31470, 31525, 31580, 31633, 31684, 31735,

    31785, 31833, 31880, 31926, 31970, 32014, 32056, 32097,

    32137, 32176, 32213, 32249, 32284, 32318, 32350, 32382,

    32412, 32441, 32468, 32495, 32520, 32544, 32567, 32588,

    32609, 32628, 32646, 32662, 32678, 32692, 32705, 32717,

    32727, 32736, 32744, 32751, 32757, 32761, 32764, 32766,

    32767, 32766, 32764, 32761, 32757, 32751, 32744, 32736,

    32727, 32717, 32705, 32692, 32678, 32662, 32646, 32628,

    32609, 32588, 32567, 32544, 32520, 32495, 32468, 32441,

    32412, 32382, 32350, 32318, 32284, 32249, 32213, 32176,

    32137, 32097, 32056, 32014, 31970, 31926, 31880, 31833,

    31785, 31735, 31684, 31633, 31580, 31525, 31470, 31413,

    31356, 31297, 31236, 31175, 31113, 31049, 30984, 30918,

    30851, 30783, 30713, 30643, 30571, 30498, 30424, 30349,

    30272, 30195, 30116, 30036, 29955, 29873, 29790, 29706,

    29621, 29534, 29446, 29358, 29268, 29177, 29085, 28992,

    28897, 28802, 28706, 28608, 28510, 28410, 28309, 28208,

    28105, 28001, 27896, 27790, 27683, 27575, 27466, 27355,

    27244, 27132, 27019, 26905, 26789, 26673, 26556, 26437,

    26318, 26198, 26077, 25954, 25831, 25707, 25582, 25456,

    25329, 25201, 25072, 24942, 24811, 24679, 24546, 24413,

    24278, 24143, 24006, 23869, 23731, 23592, 23452, 23311,

    23169, 23027, 22883, 22739, 22594, 22448, 22301, 22153,

    22004, 21855, 21705, 21554, 21402, 21249, 21096, 20942,

    20787, 20631, 20474, 20317, 20159, 20000, 19840, 19680,

    19519, 19357, 19194, 19031, 18867, 18702, 18537, 18371,

    18204, 18036, 17868, 17699, 17530, 17360, 17189, 17017,

    16845, 16672, 16499, 16325, 16150, 15975, 15799, 15623,

    15446, 15268, 15090, 14911, 14732, 14552, 14372, 14191,

    14009, 13827, 13645, 13462, 13278, 13094, 12909, 12724,

    12539, 12353, 12166, 11980, 11792, 11604, 11416, 11227,

    11038, 10849, 10659, 10469, 10278, 10087, 9895, 9703,

    9511, 9319, 9126, 8932, 8739, 8545, 8351, 8156,

    7961, 7766, 7571, 7375, 7179, 6982, 6786, 6589,

    6392, 6195, 5997, 5799, 5601, 5403, 5205, 5006,

    4807, 4608, 4409, 4210, 4011, 3811, 3611, 3411,

    3211, 3011, 2811, 2610, 2410, 2209, 2009, 1808,

    1607, 1406, 1206, 1005, 804, 603, 402, 201,

    0, -201, -402, -603, -804, -1005, -1206, -1406,

    -1607, -1808, -2009, -2209, -2410, -2610, -2811, -3011,

    -3211, -3411, -3611, -3811, -4011, -4210, -4409, -4608,

    -4807, -5006, -5205, -5403, -5601, -5799, -5997, -6195,

    -6392, -6589, -6786, -6982, -7179, -7375, -7571, -7766,

    -7961, -8156, -8351, -8545, -8739, -8932, -9126, -9319,

    -9511, -9703, -9895, -10087, -10278, -10469, -10659, -10849,

    -11038, -11227, -11416, -11604, -11792, -11980, -12166, -12353,

    -12539, -12724, -12909, -13094, -13278, -13462, -13645, -13827,

    -14009, -14191, -14372, -14552, -14732, -14911, -15090, -15268,

    -15446, -15623, -15799, -15975, -16150, -16325, -16499, -16672,

    -16845, -17017, -17189, -17360, -17530, -17699, -17868, -18036,

    -18204, -18371, -18537, -18702, -18867, -19031, -19194, -19357,

    -19519, -19680, -19840, -20000, -20159, -20317, -20474, -20631,

    -20787, -20942, -21096, -21249, -21402, -21554, -21705, -21855,

    -22004, -22153, -22301, -22448, -22594, -22739, -22883, -23027,

    -23169, -23311, -23452, -23592, -23731, -23869, -24006, -24143,

    -24278, -24413, -24546, -24679, -24811, -24942, -25072, -25201,

    -25329, -25456, -25582, -25707, -25831, -25954, -26077, -26198,

    -26318, -26437, -26556, -26673, -26789, -26905, -27019, -27132,

    -27244, -27355, -27466, -27575, -27683, -27790, -27896, -28001,

    -28105, -28208, -28309, -28410, -28510, -28608, -28706, -28802,

    -28897, -28992, -29085, -29177, -29268, -29358, -29446, -29534,

    -29621, -29706, -29790, -29873, -29955, -30036, -30116, -30195,

    -30272, -30349, -30424, -30498, -30571, -30643, -30713, -30783,

    -30851, -30918, -30984, -31049, -31113, -31175, -31236, -31297,

    -31356, -31413, -31470, -31525, -31580, -31633, -31684, -31735,

    -31785, -31833, -31880, -31926, -31970, -32014, -32056, -32097,

    -32137, -32176, -32213, -32249, -32284, -32318, -32350, -32382,

    -32412, -32441, -32468, -32495, -32520, -32544, -32567, -32588,

    -32609, -32628, -32646, -32662, -32678, -32692, -32705, -32717,

    -32727, -32736, -32744, -32751, -32757, -32761, -32764, -32766,

};*/

const myint2 Sinewave[N_WAVE-N_WAVE/4] = {
0,13176700,26352905,39528118,52701842,65873582,79042842,92209126,
105371939,118530784,131685167,144834592,157978564,171116588,184248170,197372814,
210490028,223599317,236700188,249792147,262874702,275947359,289009628,302061015,
315101030,328129181,341144979,354147933,367137553,380113351,393074838,406021525,
418952927,431868555,444767923,457650546,470515939,483363617,496193097,509003896,
521795531,534567520,547319384,560050641,572760813,585449421,598115986,610760034,
623381086,635978668,648552307,661101527,673625858,686124827,698597964,711044799,
723464864,735857691,748222813,760559765,772868082,785147302,797396961,809616599,
821805755,833963970,846090788,858185750,870248403,882278291,894274962,906237964,
918166847,930061161,941920459,953744295,965532222,977283798,988998580,1000676127,
1012315998,1023917757,1035480966,1047005189,1058489994,1069934947,1081339618,1092703577,
1104026396,1115307649,1126546912,1137743761,1148897775,1160008533,1171075618,1182098613,
1193077102,1204010672,1214898913,1225741413,1236537765,1247287562,1257990399,1268645874,
1279253585,1289813133,1300324121,1310786152,1321198832,1331561771,1341874577,1352136862,
1362348240,1372508327,1382616739,1392673097,1402677021,1412628136,1422526066,1432370439,
1442160884,1451897033,1461578518,1471204976,1480776044,1490291362,1499750572,1509153316,
1518499242,1527787998,1537019233,1546192600,1555307754,1564364352,1573362052,1582300517,
1591179409,1599998394,1608757139,1617455317,1626092598,1634668657,1643183173,1651635823,
1660026291,1668354260,1676619416,1684821448,1692960048,1701034909,1709045728,1716992201,
1724874032,1732690921,1740442576,1748128704,1755749017,1763303226,1770791049,1778212202,
1785566406,1792853385,1800072865,1807224572,1814308239,1821323599,1828270386,1835148340,
1841957203,1848696716,1855366628,1861966686,1868496642,1874956250,1881345267,1887663453,
1893910570,1900086382,1906190657,1912223165,1918183680,1924071975,1929887831,1935631028,
1941301349,1946898581,1952422514,1957872940,1963249653,1968552451,1973781133,1978935505,
1984015371,1989020539,1993950823,1998806035,2003585994,2008290518,2012919432,2017472561,
2021949733,2026350781,2030675537,2034923839,2039095528,2043190447,2047208440,2051149358,
2055013051,2058799374,2062508184,2066139342,2069692712,2073168159,2076565552,2079884764,
2083125670,2086288147,2089372077,2092377343,2095303833,2098151435,2100920044,2103609554,
2106219865,2108750878,2111202498,2113574632,2115867191,2118080090,2120213244,2122266573,
2124240001,2126133452,2127946855,2129680143,2131333249,2132906113,2134398673,2135810875,
2137142665,2138393993,2139564811,2140655077,2141664748,2142593786,2143442158,2144209830,
2144896774,2145502964,2146028377,2146472994,2146836797,2147119773,2147321911,2147443204,
2147483647,2147443239,2147321981,2147119878,2146836937,2146473168,2146028587,2145503209,
2144897054,2144210145,2143442507,2142594171,2141665167,2140655531,2139565300,2138394516,
2137143223,2135811468,2134399301,2132906775,2131333947,2129680875,2127947622,2126134253,
2124240837,2122267444,2120214149,2118081030,2115868166,2113575641,2111203541,2108751955,
2106220977,2103610700,2100921224,2098152650,2095305081,2092378626,2089373393,2086289498,
2083127054,2079886183,2076567004,2073169645,2069694232,2066140896,2062509771,2058800994,
2055014705,2051151045,2047210161,2043192201,2039097316,2034925660,2030677391,2026352667,
2021951653,2017474514,2012921418,2008292537,2003588044,1998808118,1993952939,1989022688,
1984017551,1978937718,1973783379,1968554728,1963251962,1957875281,1952424887,1946900986,
1941303785,1935633496,1929890330,1924074506,1918186242,1912225759,1906193282,1900089038,
1893913256,1887666171,1881348015,1874959028,1868499450,1861969525,1855369497,1848699616,
1841960132,1835151300,1828273376,1821326618,1814311288,1807227651,1800075972,1792856522,
1785569572,1778215397,1770794273,1763306479,1755752298,1748132014,1740445914,1732694288,
1724877426,1716995624,1709049178,1701038388,1692963554,1684824982,1676622977,1668357848,
1660029906,1651639466,1643186842,1634672353,1626096320,1617459065,1608760914,1600002194,
1591183236,1582304370,1573365931,1564368256,1555311684,1546196555,1537023212,1527792002,
1518503272,1509157370,1499754650,1490295465,1480780172,1471209127,1461582693,1451901231,
1442165106,1432374685,1422530335,1412632428,1402681336,1392677435,1382621099,1372512709,
1362352645,1352141289,1341879026,1331566242,1321203325,1310790666,1300328656,1289817690,
1279258162,1268650472,1257995018,1247292201,1236542424,1225746092,1214903612,1204015391,
1193081840,1182103370,1171080395,1160013329,1148902589,1137748594,1126551764,1115312519,
1104031284,1092708482,1081344541,1069939888,1058494952,1047010165,1035485958,1023922766,
1012321024,1000681169,989003638,977288872,965537312,953749400,941925580,930066298,
918171998,906243130,894280143,882283486,870253613,858190974,846096025,833969222,
821811020,809621877,797402252,785152606,772873399,760565094,748228154,735863044,
723470229,711050176,698603353,686130227,673631269,661106949,648557739,635984111,
623386539,610765497,598121459,585454903,572766305,560056143,547324894,534573040,
521801059,509009432,496198642,483369170,470521499,457656114,444773498,431874137,
418958516,406027121,393080440,380118959,367143168,354153553,341150605,328134813,
315106667,302066657,289015274,275953011,262880358,249797807,236705852,223604985,
210495699,197378489,184253847,171122268,157984247,144840277,131690855,118536474,
105377631,92214820,79048537,65879278,52707539,39533815,26358603,13182399,
5698,-13171002,-26347207,-39522420,-52696145,-65867886,-79037147,-92203433,
-105366247,-118525094,-131679479,-144828906,-157972881,-171110907,-184242492,-197367140,
-210484357,-223593650,-236694524,-249786487,-262869046,-275941708,-289003981,-302055373,
-315095393,-328123550,-341139353,-354142312,-367131938,-380107742,-393069235,-406015930,
-418947338,-431862972,-444762348,-457644979,-470510379,-483358065,-496187553,-508998360,
-521790003,-534562001,-547313874,-560045140,-572755321,-585443938,-598110513,-610754570,
-623375633,-635973226,-648546874,-661096106,-673620447,-686119427,-698592576,-711039422,
-723459499,-735852337,-748217471,-760554436,-772862766,-785141998,-797391670,-809611321,
-821800490,-833958719,-846085550,-858180527,-870243193,-882273096,-894269781,-906232798,
-918161695,-930056025,-941915338,-953739189,-965527132,-977278724,-988993522,-1000671084,
-1012310973,-1023912748,-1035475973,-1047000214,-1058485036,-1069930006,-1081334694,-1092698671,
-1104021508,-1115302780,-1126542061,-1137738928,-1148892960,-1160003738,-1171070841,-1182093855,
-1193072364,-1204005954,-1214894214,-1225736734,-1236533106,-1247282923,-1257985781,-1268641276,
-1279249008,-1289808577,-1300319586,-1310781638,-1321194340,-1331557300,-1341870128,-1352132435,
-1362343835,-1372503944,-1382612379,-1392668759,-1402672706,-1412623844,-1422521797,-1432366193,
-1442156662,-1451892834,-1461574343,-1471200825,-1480771917,-1490287259,-1499746493,-1509149262,
-1518495213,-1527783993,-1537015253,-1546188645,-1555303825,-1564360448,-1573358174,-1582296664,
-1591175582,-1599994593,-1608753365,-1617451568,-1626088876,-1634664962,-1643179504,-1651632181,
-1660022676,-1668350672,-1676615855,-1684817915,-1692956542,-1701031431,-1709042277,-1716988779,
-1724870637,-1732687555,-1740439238,-1748125395,-1755745736,-1763299974,-1770787825,-1778209007,
-1785563240,-1792850249,-1800069757,-1807221494,-1814305191,-1821320579,-1828267397,-1835145381,
-1841954273,-1848693817,-1855363758,-1861963846,-1868493833,-1874953472,-1881342520,-1887660736,
-1893907884,-1900083727,-1906188033,-1912220572,-1918181117,-1924069444,-1929885331,-1935628560,
-1941298912,-1946896177,-1952420141,-1957870599,-1963247343,-1968550173,-1973778888,-1978933292,
-1984013190,-1989018391,-1993948707,-1998803952,-2003583943,-2008288500,-2012917447,-2017470609,
-2021947814,-2026348894,-2030673683,-2034922019,-2039093741,-2043188692,-2047206719,-2051147670,
-2055011396,-2058797753,-2062506597,-2066137789,-2069691192,-2073166672,-2076564100,-2079883345,
-2083124285,-2086286796,-2089370760,-2092376060,-2095302584,-2098150221,-2100918864,-2103608408,
-2106218753,-2108749800,-2111201455,-2113573623,-2115866217,-2118079150,-2120212339,-2122265703,
-2124239165,-2126132650,-2127946088,-2129679411,-2131332552,-2132905450,-2134398045,-2135810282,
-2137142106,-2138393469,-2139564322,-2140654623,-2141664329,-2142593402,-2143441808,-2144209515,
-2144896494,-2145502719,-2146028167,-2146472819,-2146836657,-2147119668,-2147321841,-2147443169
};
/*

FIX_MPY() - fixed-point multiplication & scaling.

Substitute inline assembly for hardware-specific

optimization suited to a particluar DSP processor.

Scaling ensures that result remains 16-bit.

*/

inline short FIX_MPY(short a, short b)

{

    /* shift right one less bit (i.e. 15-1) */
//ligb 20170417把下面int c改为unsigned char c，以免在有些机器符号不正确。决定放弃修改，因为。
    int c = ((int)a * (int)b) >> 14;

    /* last bit shifted out = rounding-bit */

    b = c & 0x01;

    /* last shift + rounding bit */

    a = (c >> 1) + b;

    return a;

}

/*

fix_fft() - perform forward/inverse fast Fourier transform.

fr[n],fi[n] are real and imaginary arrays, both INPUT AND

RESULT (in-place FFT), with 0 <= n < 2**m; set inverse to

0 for forward transform (FFT), or 1 for iFFT.

*/

int fix_fft(myint fr[], myint fi[], short m, short inverse)

{
    myint nHalfInt;//用于判断移位，它是8位或16位有符号数最大值的一半 ligb 20170417
    if (sizeof(myint)==1 )
        nHalfInt=127;
    if (sizeof(myint)==2)
        nHalfInt=16383;

    int mr, nn, i, j,j1, l,  istep, n, scale, shift;
    unsigned char k;
    unsigned char k1;//k1  is different power of 2 for N_WAVE to current FFT points number
//下句原为k=LOG2_N_WAVE-1 ;//为了适应1024以下点数FFT而改写成依赖FFT点数n的变量
    k=LOG2_N_WAVE-1 ;
    //k = log(n)/log(2)-1;//ligb 20170417
    k1=LOG2_N_WAVE-k-2;//k1是j在旋转因子中索引应该扩大的2的幂指数

    short qr, qi, tr, ti, wr, wi;

    n = 1 << m;

    /* max FFT size = N_WAVE */

    if (n > N_WAVE)

        return -1;

    mr = 0;

    nn = n - 1;

    scale = 0;

    /* 抽取decimation in time - re-order data */

    for (m=1; m<=nn; ++m)
    {

        l = n;

        do
        {

            l >>= 1;

        }
        while (mr+l > nn);

        mr = (mr & (l-1)) + l;

        if (mr <= m)

            continue;

        tr = fr[m];

        fr[m] = fr[mr];

        fr[mr] = tr;

        ti = fi[m];

        fi[m] = fi[mr];

        fi[mr] = ti;

    }

    l = 1;

    while (l < n)
    {

        if (inverse)
        {

            /* variable scaling, depending upon data */

            shift = 0;

            for (i=0; i<n; ++i)
            {

                j = fr[i];

                if (j < 0)

                    j = -j;

                m = fi[i];

                if (m < 0)

                    m = -m;

                if (j > nHalfInt || m > nHalfInt)
                {

                    shift = 1;

                    break;

                }

            }

            if (shift)

                ++scale;

        }
        else
        {

            /*

            fixed scaling, for proper normalization --

            there will be log2(n) passes, so this results

            in an overall factor of 1/n, distributed to

            maximize arithmetic accuracy.

            */

            shift = 1;

        }

        /*

        it may not be obvious, but the shift will be

        performed on each data point exactly once,

        during this pass.

        */

        istep = l << 1;

        for (m=0; m<l; ++m)
        {

            j = m << k;
            j1=j<<k1;//ligb 20170419

            /* 0 <= j < N_WAVE/2 */

            wr = Sinewave[j1+N_WAVE/4];

            wi = -Sinewave[j1];

            if (inverse)

                wi = -wi;

            if (shift)
            {

                wr >>= 1;

                wi >>= 1;

            }

            for (i=m; i<n; i+=istep)
            {

                j = i + l;

                tr = FIX_MPY(wr,fr[j]) - FIX_MPY(wi,fi[j]);

                ti = FIX_MPY(wr,fi[j]) + FIX_MPY(wi,fr[j]);

                qr = fr[i];

                qi = fi[i];

                if (shift)
                {

                    qr >>= 1;

                    qi >>= 1;

                }

                fr[j] = qr - tr;

                fi[j] = qi - ti;

                fr[i] = qr + tr;

                fi[i] = qi + ti;

            }

        }

        --k;

        l = istep;

    }

    return scale;

}

/*

fix_fftr() - forward/inverse FFT on array of real numbers.

Real FFT/iFFT using half-size complex FFT by distributing

even/odd samples into real/imaginary arrays respectively.

In order to save data space (i.e. to avoid two arrays, one

for real, one for imaginary samples), we proceed in the

following two steps: a) samples are rearranged in the real

array so that all even samples are in places 0-(N/2-1) and

all imaginary samples in places (N/2)-(N-1), and b) fix_fft

is called with fr and fi pointing to index 0 and index N/2

respectively in the original array. The above guarantees

that fix_fft "sees" consecutive real samples as alternating

real and imaginary samples in the complex array.

*/

int fix_fftr(myint f[], short m, short inverse)

{

    int i, scale = 0;
    int N = 1<<(m-1);

    short tt, *fr=f, *fi=&f[N];

    if (inverse)

        scale = fix_fft(fr, fi, m-1, inverse);

    for (i=1; i<N; i+=2)
    {

        tt = f[N+i-1];

        f[N+i-1] = f[i];

        f[i] = tt;

    }

    if (! inverse)

        scale = fix_fft(fi, fr, m-1, inverse);

    return scale;

}
